/*
 *  Copyright (c) 2015-2018 Jose Osvaldo Suarez Domingos
 *
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
 * 
 */

///////////////////////////////////////////////////////////////////////////////

//Sustitucion para std::enable_if
template<bool B, class T = void>
struct enable_if {};
 
template<class T>
struct enable_if<true, T> {typedef T type;};

///////////////////////////////////////////////////////////////////////////////

template<class T, T v>
struct integral_constant
{
    static constexpr T value = v;
    typedef T value_type;
    typedef integral_constant type;
    constexpr operator value_type() const noexcept {return value;}
    constexpr value_type operator()() const noexcept {return value;}
};

///////////////////////////////////////////////////////////////////////////////

typedef integral_constant<bool, true> true_type;
typedef integral_constant<bool, false> false_type;
  
///////////////////////////////////////////////////////////////////////////////

//Sustitucion para std::is_same
template<class T, class U>
struct is_same : false_type {};
 
template<class T>
struct is_same<T, T> : true_type {};

///////////////////////////////////////////////////////////////////////////////

//Clase que permite saber si un tipo es "big_uint".
//template<typename T>
//struct is_big_uint : false_type {};

//template<unsigned int A, typename B>
//struct is_big_uint<big_uint<A, B>> : true_type {};

template<typename T>
struct is_big_uint
{
    template <typename C>
    static constexpr bool test(decltype(C::make_big_uint) *)
    {
        return true;
        //return C::is_trivially_copyable;
    }

    template <typename C>
    static constexpr bool test(...)
    {
        return false;
    }

    static constexpr bool value = test<T>(nullptr);
};

///////////////////////////////////////////////////////////////////////////////

template<typename From,
         typename To>
struct conversion_loses_precision : true_type {};

template<typename T>
struct conversion_loses_precision<T, T> : false_type {};

template<> struct conversion_loses_precision<uint8_t, uint64_t> : false_type {};
template<> struct conversion_loses_precision<uint16_t, uint64_t> : false_type {};
template<> struct conversion_loses_precision<uint32_t, uint64_t> : false_type {};
template<> struct conversion_loses_precision<uint8_t, uint32_t> : false_type {};
template<> struct conversion_loses_precision<uint16_t, uint32_t> : false_type {};
template<> struct conversion_loses_precision<uint8_t, uint16_t> : false_type {};

template<> struct conversion_loses_precision<int8_t, uint64_t> : false_type {};
template<> struct conversion_loses_precision<int16_t, uint64_t> : false_type {};
template<> struct conversion_loses_precision<int32_t, uint64_t> : false_type {};
template<> struct conversion_loses_precision<int8_t, uint32_t> : false_type {};
template<> struct conversion_loses_precision<int16_t, uint32_t> : false_type {};
template<> struct conversion_loses_precision<int8_t, uint16_t> : false_type {};
template<> struct conversion_loses_precision<int8_t, uint8_t> : false_type {};

///////////////////////////////////////////////////////////////////////////////

template<typename T>
struct is_negative_helper
{
    static constexpr bool is_negative(const T & n)
    {
        return n < static_cast<T>(0);
    }
};
template<> struct is_negative_helper<uint8_t>
{static constexpr bool is_negative(const uint8_t &) {return false;}};
template<> struct is_negative_helper<uint16_t>
{static constexpr bool is_negative(const uint16_t &) {return false;}};
template<> struct is_negative_helper<uint32_t>
{static constexpr bool is_negative(const uint32_t &) {return false;}};
template<> struct is_negative_helper<uint64_t>
{static constexpr bool is_negative(const uint64_t &) {return false;}};

template<typename T>
struct is_unsigned : false_type {};
template<> struct is_unsigned<uint8_t> : false_type {};
template<> struct is_unsigned<uint16_t> : false_type {};
template<> struct is_unsigned<uint32_t> : false_type {};
template<> struct is_unsigned<uint64_t> : false_type {};

template<typename T>
constexpr bool is_negative(const T & n)
{
    return is_negative_helper<T>::is_negative(n);
}

///////////////////////////////////////////////////////////////////////////////

//Desplazamiento a la izquierda
//Nota: La comprobacion de que no se excede el tamano del tipo de datos
//se hace en tiempo de compilacion para evitar molestas advertencias del
//compilador.
template<unsigned int Shift_Count, typename T>
static T shift_left(const T & value)
{
    constexpr bool enough_bits = (sizeof(T) * 8 > Shift_Count);
    using enough_bits_t = integral_constant<bool, enough_bits>;
    return shift_left<Shift_Count, T>(enough_bits_t{}, value);
}

template<unsigned int Shift_Count, typename T>
static T shift_left(false_type, const T &)
{
    return static_cast<T>(0);
}

template<unsigned int Shift_Count, typename T>
static T shift_left(true_type, const T & value)
{
    return value << Shift_Count;
}

///////////////////////////////////////////////////////////////////////////////

//Desplazamiento a la derecha.
//Nota: La comprobacion de que no se excede el tamano del tipo de datos
//se hace en tiempo de compilacion para evitar molestas advertencias del
//compilador.
template<unsigned int Shift_Count, typename T>
static T shift_right(const T & value)
{
    constexpr bool enough_bits = (sizeof(T) * 8 > Shift_Count);
    using enough_bits_t = integral_constant<bool, enough_bits>;
    return shift_right<Shift_Count, T>(enough_bits_t{}, value);
}

template<unsigned int Shift_Count, typename T>
static T shift_right(false_type, const T &)
{
    return static_cast<T>(0);
}

template<unsigned int Shift_Count, typename T>
static T shift_right(true_type, const T & value)
{
    return value >> Shift_Count;
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
void zero_blocks(T * result,
                 const unsigned int blocks)
{
    const T zero = static_cast<T>(0);

    for (unsigned int i = 0; i < blocks; i++)
        result[i] = zero;
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
void fill_blocks(T * result,
                 const unsigned int blocks)
{
    const T all = ~static_cast<T>(0);

    for (unsigned int i = 0; i < blocks; i++)
        result[i] = all;
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
void copy_blocks(T * to,
                 const T * from,
                 const unsigned int blocks)
{
    for (unsigned int i = 0; i < blocks; i++)
        to[i] = from[i];
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
void not_blocks(T * result,
                const T * data,
                const unsigned int blocks)
{
    for (unsigned int i = 0; i < blocks; i++)
        result[i] = ~data[i];
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
void and_blocks(T * result,
                const T * a,
                const T * b,
                const unsigned int blocks)
{
    for (unsigned int i = 0; i < blocks; i++)
        result[i] = a[i] & b[i];
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
void or_blocks(T * result,
               const T * a,
               const T * b,
               const unsigned int blocks)
{
    for (unsigned int i = 0; i < blocks; i++)
        result[i] = a[i] | b[i];
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
void xor_blocks(T * result,
                const T * a,
                const T * b,
                const unsigned int blocks)
{
    for (unsigned int i = 0; i < blocks; i++)
        result[i] = a[i] ^ b[i];
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
T addc(const T a,
       const T b,
       T & carry)
{
#if 0    
    const unsigned int bits = sizeof(T) * 8;
    const T one = static_cast<T>(1);
    const T mask_high = one << (bits - 1);
    const T mask_low = mask_high - 1;

    const T a_low = a & mask_low;
    const T b_low = b & mask_low;
    const T r_low = a_low + b_low + carry;
    const T r_high = r_low >> (bits - 1);
    const T a_high = a >> (bits - 1);
    const T b_high = b >> (bits - 1);
    const T high = r_high + a_high + b_high;

    carry = high >= 2 ? 1 : 0;

    const T r = (r_low & mask_low) +
                ((high & 0x1) << (bits - 1));
                constexpr T zero = static_cast<T>(0);
#else
    T r = a + b;
    const bool overflow = r < a || r < b;
    if (carry)
        r++;
    carry = overflow || (carry && r == 0);
#endif

    return r;
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
void addc_blocks(T * result,
                 const T * a,
                 const T * b,
                 T & carry,
                 const unsigned int blocks)
{
    for (unsigned int i = 0; i < blocks; i++)
        result[i] = addc<T>(a[i], b[i], carry);
}

///////////////////////////////////////////////////////////////////////////////

//Propaga el acarreo.
template<typename T>
inline
void inc_blocks(T * result,
                T & carry,
                const unsigned int blocks)
{
    for (unsigned int i = 0; i < blocks && carry; i++)
        result[i] = addc<T>(result[i], static_cast<T>(0), carry);
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
void neg_blocks(T * result,
                const T * data,
                T & carry,
                const unsigned int blocks)
{
    //not_blocks
    //inc_blocks

    for (unsigned int i = 0; i < blocks; i++)
        result[i] = addc<T>(~data[i], static_cast<T>(0), carry);
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
void subc_blocks(T * result,
                 const T * a,
                 const T * b,
                 T & carry,
                 const unsigned int blocks)
{
    for (unsigned int i = 0; i < blocks; i++)
        result[i] = addc<T>(a[i], ~b[i], carry);
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
void dec_blocks(T * result,
                T & carry,
                const unsigned int blocks)
{
    for (unsigned int i = 0; i < blocks && !carry; i++)
        result[i] = addc<T>(result[i], ~static_cast<T>(0), carry);
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
void shift_left_big(T * result,
                    const T * data,
                    const unsigned int n,
                    const unsigned int blocks)
{
    const T zero = static_cast<T>(0);
    const unsigned int bits_per_block = sizeof(T) * 8u;
    const unsigned int blocks_skip = n / bits_per_block;
    const unsigned int bits_skip = n % bits_per_block;

    if (blocks_skip < blocks)
        zero_blocks(result, blocks_skip);
    else
        zero_blocks(result, blocks);

    if (bits_skip == 0)
        for (unsigned int i = blocks_skip; i < blocks; i++)
            result[i] = data[i - blocks_skip];
    else
    {
        T result_i = zero;
        T data_j;
        T data_j_1 = zero;
        for (unsigned int i = blocks_skip; i < blocks; i++)
        {
            data_j = data_j_1;
            result_i = data_j >> (bits_per_block - bits_skip);
            data_j_1 = data[i - blocks_skip];
            result_i |= data_j_1 << bits_skip;
            result[i] = result_i;
        }
    }
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
void shift_right_big(T * result,
                     const T * data,
                     const unsigned int n,
                     const unsigned int blocks)
{
    //const T zero = static_cast<T>(0);
    const unsigned int bits_per_block = sizeof(T) * 8u;
    const unsigned int blocks_skip = n / bits_per_block;
    const unsigned int bits_skip = n % bits_per_block;

    if (blocks_skip >= blocks)
    {
        zero_blocks(result, blocks);
        return;
    }

    if (bits_skip == 0)
        for (unsigned int i = 0; i < blocks - blocks_skip; i++)
            result[i] = data[i + blocks_skip];
    else
    {
        for (unsigned int i = 0; i < blocks - blocks_skip; i++)
        {
            result[i] = data[i + blocks_skip] >> bits_skip;

            if (i + blocks_skip + 1 < blocks)
                result[i] |= data[i + blocks_skip + 1] <<
                             (bits_per_block - bits_skip);
        }
    }

    zero_blocks(result + blocks - blocks_skip,
                blocks_skip);
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
bool less(const T * a,
          const T * b,
          const unsigned int blocks)
{
    unsigned int i = blocks;
    while(i > 0)
    {
        const unsigned int j = i - 1;

        if (a[j] < b[j])
            return true;
        else if (b[j] < a[j])
            return false;

        i--;
    }

    return false;
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
bool equal(const T * a,
           const T * b,
           const unsigned int blocks)
{
    unsigned int i = 0;
    while(i < blocks)
    {
        if (!(a[i] == b[i]))
            return false;
        i++;
    }

    return true;
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
bool less_equal(const T * a,
                const T * b,
                const unsigned int blocks)
{
    unsigned int i = blocks;
    while(i > 0)
    {
        const unsigned int j = i - 1;

        if (a[j] < b[j])
            return true;
        else if (b[j] < a[j])
            return false;

        i--;
    }

    return true;
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline
T mul(T a,
      T b,
      T & high)
{
    constexpr unsigned int bits = sizeof(T) * 8;
    constexpr unsigned int half_bits = bits / 2;
    constexpr T mask_low = (static_cast<T>(1) << half_bits) - 1;
    constexpr T mask_high = ~mask_low;
    constexpr T zero = static_cast<T>(0);
    //constexpr T one = static_cast<T>(1);

    T a_low = a & mask_low;
    T a_high = (a & mask_high) >> half_bits;
    T b_low = b & mask_low;
    T b_high = (b & mask_high) >> half_bits;

    T low = zero;
    high = zero;

    T carry = zero;
    low = addc<T>(low,
                  a_low * b_low,
                  carry);
    high += carry;

    carry = zero;
    low = addc<T>(low,
                  (a_high * b_low) << half_bits,
                  carry);
    high += carry;

    carry = zero;
    low = addc<T>(low,
                  (a_low * b_high) << half_bits,
                  carry);
    high += carry;

    high += (a_high * b_low) >> half_bits;
    high += (a_low * b_high) >> half_bits;
    high += a_high * b_high;

    return low;
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
void fma2_kernel(T * data,
                 T a,
                 const T * b_data,
                 const unsigned int blocks_data,
                 const unsigned int blocks_b)
{
    using base_type = T;
    const base_type zero = static_cast<base_type>(0);
    const base_type all = ~zero;

    base_type tmp0;
    base_type tmp1 = data[0];
    unsigned int j;
    for (j = 0; j < blocks_b; j++)
    {
        unsigned int k = j;
        if (k == blocks_data) break;
                 
        tmp0 = tmp1;

        base_type high = zero;
        base_type low = mul<T>(a, b_data[k], high);

        T carry = zero;
        data[k] = addc<T>(tmp0, low, carry);
        k++;
        if (k == blocks_data) break;
        tmp1 = addc<T>(data[k], high, carry);
        k++;

        if (carry == zero) continue;

        while(k < blocks_data &&
              data[k] == all)
        {
            data[k]++;
            k++;
        }
        if (k < blocks_data) data[k]++;
    }//end for j

    if (j == blocks_b &&
        j < blocks_data)
        data[j] = tmp1;
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
void fma2(T * data,
          const unsigned int blocks_data,
          const T * a_data,
          const unsigned int blocks_a,
          const T * b_data,
          const unsigned int blocks_b)
{
    using base_type = T;

    for (unsigned int i = 0;
         i < blocks_a && i < blocks_data;
         i++)
    {
        const base_type a = a_data[i];

        fma2_kernel<base_type>(data + i,
                               a,
                               b_data,
                               blocks_data - i,
                               blocks_b);
    }//end for i
}

///////////////////////////////////////////////////////////////////////////////

template<typename Base_Type,
         unsigned int N,
         bool is_pow2 = N == 1 ||
                        N == 2 ||
                        N == 4 ||
                        N == 8 ||
                        N == 16 ||
                        N == 32>
struct addmul_1;

template<typename Base_Type,
         unsigned int N>
struct addmul_1<Base_Type, N, false>
{
    using base_type = Base_Type;

    static void operate(base_type * data,
                        const base_type a0,
                        const base_type * b,
                        base_type & tmp0)
    {
        const unsigned int burst = N >= 16 ? 16 :
                                   N >= 8  ? 8  :
                                   N >= 4  ? 4  :
                                   N >= 2  ? 2  :
                                             1;
        addmul_1<base_type, burst>::operate(data + 0,
                                            a0,
                                            b + 0,
                                            tmp0);
        addmul_1<base_type, N - burst>::operate(data + burst,
                                                a0,
                                                b + burst,
                                                tmp0);
    }
};

template<typename Base_Type,
         unsigned int N>
struct addmul_1<Base_Type, N, true>
{
    using base_type = Base_Type;

    static void operate(base_type * data,
                        const base_type a0,
                        const base_type * b,
                        base_type & tmp0)
    {
        const unsigned int N2 = N / 2;
        for (unsigned int i = 0; i < 2; i++)
            addmul_1<base_type, N2>::operate(data + i * N2,
                                             a0,
                                             b + i * N2,
                                             tmp0);
    }
};

template<typename Base_Type>
struct addmul_1<Base_Type, 1>
{
    using base_type = Base_Type;
    static const base_type zero = static_cast<base_type>(0);

    static void operate(base_type * data,
                        const base_type a0,
                        const base_type * b,
                        base_type & tmp0)
    {
        base_type low0,
                  high0,
                  carry, b0;
        b0 = b[0];
        low0 = mul<base_type>(a0, b0, high0);
        carry = zero;
        tmp0 = addc<base_type>(tmp0, low0, carry);
        high0 += carry;
        carry = zero;
        data[0] = addc<base_type>(data[0], tmp0, carry);
        //tmp0 = addc<base_type>(tmp1, high0, carry);
        //tmp1 = carry;
        tmp0 = high0 + carry;

        //Se puede demostrar que este ultimo acarreo siempre es 0.
        //
        //Partimos de que la primera vez (tmp1 == 0),
        //por tanto, para que se genere el primer acarreo
        //(tmp1) no cuenta.
        //
        //La unica posibilidad que queda es que
        //(high0 == 0xFFFF...) y (carry == 1) simultaneamente.
        //Se va a ver que esto no es posible.
        //
        //Concretamente,
        //se va a demostrar que si (high0 == 0xFFFF), (tmp == 0),
        //haciendo que la suma de (data[0] + tmp0) nunca genere
        //acarreo, pues no hay que olvidar que el acarreo de
        //entrada de esa suma es 0.
        //
        //Se comprueba que el valor maximo del producto de dos
        //enteros de N bits (por ejemplo, 16) es
        //(0xFFFF * 0xFFFF == 0xFFFE0001). Asi que el valor
        //maximo de la parte alta es (0xFFFE).
        //
        //Por tanto la operacion "high0 += carry" puede alcanzar
        //el valor (0xFFFF), pues (carry == 0 || carry == 1).
        //
        //Sin embargo, para que (high0 == 0xFFFE + carry) es
        //necesario que (low0 == 0x0001). Podriamos tener valores
        //distintos para (low0), mayores, pero eso implica que
        //(high0 < 0xFFFE + carry).
        //
        //Acaso es importante que (low0 == 0x0001)?
        //
        //Si, pues la unica forma de que el acarreo sumado a
        //(high0) sea (carry == 1) es que la suma (tmp0 + low0)
        //genere acarreo, y si (low0 == 0x0001) es necesario
        //que (tmp0 == 0xFFFF), por tanto la suma sera
        //(0xFFFF + 0x0001 == 0x0000) y genera acarreo.
        //
        //Es decir, (low0) podria ser mayor que (0x0001), pero
        //solo cuando (high0 < 0xFFFE), asi que cuando
        //(tmp0 + low0) genera acarreo resulta que o bien
        //(high0 == 0xFFFE + 1 && tmp0 == 0x0000) o bien
        //(high0 < 0xFFFE + 1).
    }
};

///////////////////////////////////////////////////////////////////////////////

template<typename Base_Type,
         unsigned int N,
         bool is_pow2 = N == 1 ||
                        N == 2 ||
                        N == 4 ||
                        N == 8 ||
                        N == 16 ||
                        N == 32>
struct addmul_2;

template<typename Base_Type,
         unsigned int N>
struct addmul_2<Base_Type, N, false>
{
    using base_type = Base_Type;

    static void operate(base_type * data,
                        const base_type a0,
                        const base_type a1,
                        const base_type * b,
                        base_type & tmp0,
                        base_type & tmp1)
    {
        const unsigned int burst = N >= 16 ? 16 :
                                   N >= 8  ? 8  :
                                   N >= 4  ? 4  :
                                   N >= 2  ? 2  :
                                             1;
        addmul_2<base_type, burst>::operate(data + 0,
                                            a0, a1,
                                            b + 0,
                                            tmp0, tmp1);
        addmul_2<base_type, N - burst>::operate(data + burst,
                                                a0, a1,
                                                b + burst,
                                                tmp0, tmp1);
    }
};

template<typename Base_Type,
         unsigned int N>
struct addmul_2<Base_Type, N, true>
{
    using base_type = Base_Type;

    static void operate(base_type * data,
                        const base_type a0,
                        const base_type a1,
                        const base_type * b,
                        base_type & tmp0,
                        base_type & tmp1)
    {
        const unsigned int N2 = N / 2;
        for (unsigned int i = 0; i < 2; i++)
            addmul_2<base_type, N2>::operate(data + i * N2,
                                             a0, a1,
                                             b + i * N2,
                                             tmp0, tmp1);
    }
};

template<typename Base_Type>
struct addmul_2<Base_Type, 1>
{
    using base_type = Base_Type;
    static const base_type zero = static_cast<base_type>(0);

    static void operate(base_type * data,
                        const base_type a0,
                        const base_type a1,
                        const base_type * b,
                        base_type & tmp0,
                        base_type & tmp1)
    {
        base_type low0, low1,
                  high0, high1,
                  carry, b0;
        b0 = b[0];
        low0 = mul<base_type>(a0, b0, high0);
        low1 = mul<base_type>(a1, b0, high1);
        carry = zero;
        tmp0 = addc<base_type>(tmp0, low0, carry);
        high0 = addc<base_type>(high0, low1, carry);
        high1 += carry;
        carry = zero;
        data[0] = addc<base_type>(data[0], tmp0, carry);
        tmp0 = addc<base_type>(tmp1, high0, carry);
        tmp1 = high1 + carry;
    }
};

///////////////////////////////////////////////////////////////////////////////

#define LAST_BLOCKS(n) \
        if (i + j + n < blocks_data)\
            data[i + j + n] = addc<base_type>(data[i + j + n],\
                                              tmp##n,\
                                              carry);
#define LAST_CARRY(n) \
        if (carry && i + j + n < blocks_data)\
            inc_blocks(data + i + j + n,\
                       carry,\
                       blocks_data - (i + j + n));

///////////////////////////////////////////////////////////////////////////////

template<typename T>
void fma8(T * data,
          const unsigned int blocks_data,
          const T * a_data,
          const unsigned int blocks_a,
          const T * b_data,
          const unsigned int blocks_b)
{
    using base_type = T;
    const base_type zero = static_cast<base_type>(0);

    auto min = [](unsigned int a, unsigned int b)
    {return a < b ? a : b;};

    const unsigned int limit = min(blocks_a, blocks_data);
    for (unsigned int i = 0; i < limit; i += 2)
    {
        const base_type a0 = a_data[i + 0];
        const base_type a1 = i + 1 < blocks_a ? a_data[i + 1] : zero;

        base_type tmp0 = zero;
        base_type tmp1 = zero;

        const unsigned int d = min(blocks_b, blocks_data - i);
        base_type * r = data + i + d;
        const base_type * b = b_data + d;
        signed int j = -static_cast<signed int>(d);
        while(j < 0)
        {
            addmul_2<base_type, 1>::
                operate(r + j,
                        a0, a1,
                        b + j,
                        tmp0, tmp1);
            j++;
        }

        base_type carry = zero;
        j += d;
        LAST_BLOCKS(0)
        LAST_BLOCKS(1)
        LAST_CARRY(2)
    }
}
#undef ADDMUL_1_LOOP
#undef ADDMUL_2_LOOP
#undef ADDMUL_4_LOOP
#undef LAST_BLOCKS

///////////////////////////////////////////////////////////////////////////////

template<typename Base_Type>
void karatsuba(Base_Type * r_data,
               const unsigned int blocks_data,
               const Base_Type * a_data,
               const Base_Type * b_data,
               const unsigned int blocks);

template<typename Base_Type>
void requete_fma(Base_Type * r,
                 const unsigned int blocks_r,
                 const Base_Type * a,
                 const unsigned int blocks_a,
                 const Base_Type * b,
                 const unsigned int blocks_b)
{
    using base_type = Base_Type;
    const unsigned int bits_per_block = sizeof(base_type) * 8u;

    auto min = [](unsigned int a, unsigned int b)
    {return a < b ? a : b;};
    const unsigned int min_blocks = min(blocks_a, blocks_b);
    if (min_blocks <= 3072 / bits_per_block)
    {
        if (blocks_a < blocks_b)
            fma8<base_type>(r, blocks_r,
                            a, blocks_a,
                            b, blocks_b);
        else
            fma8<base_type>(r, blocks_r,
                            b, blocks_b,
                            a, blocks_a);
    }
    else if (blocks_a < blocks_b)
    {
        requete_fma<base_type>(r, blocks_r,
                               b, blocks_b,
                               a, blocks_a);
    }
    else if (blocks_a > blocks_b)
    {
        requete_fma<base_type>(r, blocks_r,
                               a, blocks_b,
                               b, blocks_b);

        if (blocks_r > blocks_b)
            requete_fma<base_type>(r + blocks_b, blocks_r - blocks_b,
                                   a + blocks_b, blocks_a - blocks_b,
                                   b, blocks_b);
    }
    else if (blocks_a % 2 != 0)
    {
        const unsigned int blocks = blocks_a - 1;
        const unsigned int blocks2 = blocks * 2;

        if (blocks > 0)
            karatsuba(r, blocks_r, a, b, blocks);

        if (blocks_r > blocks)
        {
#if 0
            fma2_kernel(r + blocks,
                        a[blocks],
                        b,
                        blocks_r - blocks,
                        blocks);
            fma2_kernel(r + blocks,
                        b[blocks],
                        a,
                        blocks_r - blocks,
                        blocks);
#else
            fma2<base_type>(r + blocks, blocks_r - blocks,
                            a + blocks, 1,
                            b, blocks);
            fma2<base_type>(r + blocks, blocks_r - blocks,
                            b + blocks, 1,
                            a, blocks);
#endif
        }

        if (blocks_r > blocks2)
        {
            fma8<base_type>(r + blocks2, blocks_r - blocks2,
                            a + blocks, 1,
                            b + blocks, 1);
        }
    }
    else
    {
        karatsuba(r, blocks_r, a, b, blocks_a);
    }
}

///////////////////////////////////////////////////////////////////////////////

template<typename Base_Type,
         unsigned int A,
         unsigned int B,
         unsigned int R>
void fixed_size_multiply(Base_Type * r,
                         const Base_Type * a,
                         const Base_Type * b)
{
    using base_type = Base_Type;

    zero_blocks<base_type>(r, R);
    requete_fma<base_type>(r, R,
                           a, A,
                           b, B);
}

///////////////////////////////////////////////////////////////////////////////

template<typename Base_Type>
void variable_size_multiply(Base_Type * r,
                            const unsigned int R,
                            const Base_Type * a,
                            const unsigned int A,
                            const Base_Type * b,
                            const unsigned int B)
{
    using base_type = Base_Type;

    const unsigned int mix = A + (B << 5) + (R << 10);
    switch(mix)
    {
#define CASE(x, y, z) \
        case x + (y << 5) + (z << 10):\
            fixed_size_multiply<base_type, x, y, z>(r, a, b);\
            break;
        CASE(1, 1, 1)
        CASE(1, 1, 2)
        CASE(2, 2, 2)
        CASE(2, 2, 3)
        CASE(2, 2, 4)
        CASE(3, 3, 3)
        CASE(3, 3, 5)
        CASE(3, 3, 6)
        CASE(4, 4, 4)
        CASE(4, 4, 7)
        CASE(4, 4, 8)
        CASE(8, 8, 8)
        CASE(8, 8, 15)
        CASE(8, 8, 16)
#undef CASE
        default:
            zero_blocks<base_type>(r, R);
            requete_fma<base_type>(r, R,
                                   a, A,
                                   b, B);
    }
}

///////////////////////////////////////////////////////////////////////////////

#if 1
//PRE: blocks % 2 == 0
template<typename Base_Type>
void karatsuba(Base_Type * r_data,
               const unsigned int blocks_data,
               const Base_Type * a_data,
               const Base_Type * b_data,
               const unsigned int blocks)
{
    using base_type = Base_Type;
    const base_type zero = static_cast<base_type>(0);

    //auto max = [](unsigned int a, unsigned int b)
    //{return a > b ? a : b;};

    auto min = [](unsigned int a, unsigned int b)
    {return a < b ? a : b;};

    auto half = [](unsigned int n)
    {return n % 2 == 0 ? n / 2 : n / 2 + 1;};

    const bool high = blocks_data > blocks;
    const unsigned int factor = high ? 2 : 1;

    if (blocks > 1)
    {
        const unsigned int half_blocks = half(blocks);
        const base_type * a0 = a_data;
        const base_type * a1 = a_data + half_blocks;
        const base_type * b0 = b_data;
        const base_type * b1 = b_data + half_blocks;
/*
        printf("a0 = ");
        for (unsigned int i = half_blocks; i > 0; i--)
            printf("%.16" PRIX64 " ", a0[i - 1]);
        printf("\n");
        printf("a1 = ");
        for (unsigned int i = half_blocks; i > 0; i--)
            printf("%.16" PRIX64 " ", a1[i - 1]);
        printf("\n");
        printf("b0 = ");
        for (unsigned int i = half_blocks; i > 0; i--)
            printf("%.16" PRIX64 " ", b0[i - 1]);
        printf("\n");
        printf("b1 = ");
        for (unsigned int i = half_blocks; i > 0; i--)
            printf("%.16" PRIX64 " ", b1[i - 1]);
        printf("\n");
*/
        //base_type * c0 = r_data;
        base_type * c0 = new base_type[half_blocks * 2];
        zero_blocks<base_type>(c0, half_blocks * 2);
        requete_fma(c0, half_blocks * 2,
                    a0, half_blocks,
                    b0, half_blocks);

        base_type * c2 = new base_type[half_blocks * factor];
        zero_blocks<base_type>(c2, half_blocks * factor);
        requete_fma(c2, half_blocks * factor,
                    a1, half_blocks,
                    b1, half_blocks);

        bool sign_a = false;// = a0 < a1;
        bool sign_b = false;// = b1 < b0;
        base_type * diff_a = new base_type[blocks];
        base_type * diff_b = diff_a + half_blocks;
        base_type carry_a = 1;
        subc_blocks<base_type>(diff_a, a0, a1, carry_a, half_blocks);
        sign_a = carry_a ^ 0x1;

        base_type carry_b = 1;
        subc_blocks<base_type>(diff_b, b1, b0, carry_b, half_blocks);
        //printf("b: %u, %lu\n", sign_b, carry_b);
        sign_b = carry_b ^ 0x1;

        //printf("sign_a = %u, sign_b = %u\n", sign_a, sign_b);

        const bool sign = sign_a ^ sign_b;

        if (sign_a)
        {
            carry_a = 1;
            neg_blocks<base_type>(diff_a, diff_a, carry_a, half_blocks);
        }
        if (sign_b)
        {
            carry_b = 1;
            neg_blocks<base_type>(diff_b, diff_b, carry_b, half_blocks);
        }
/*
        printf("diff_a = ");
        for (unsigned int i = half_blocks; i > 0; i--)
            printf("%.16" PRIX64 " ", diff_a[i - 1]);
        printf("\n");
        printf("diff_b = ");
        for (unsigned int i = half_blocks; i > 0; i--)
            printf("%.16" PRIX64 " ", diff_b[i - 1]);
        printf("\n");
*/
        base_type * diff_ab = new base_type[half_blocks * factor];
        zero_blocks<base_type>(diff_ab, half_blocks * factor);
        requete_fma(diff_ab, half_blocks * factor,
                    diff_a, half_blocks,
                    diff_b, half_blocks);
/*
        printf("diff_ab = ");
        for (unsigned int i = half_blocks * factor; i > 0; i--)
            printf("%.16" PRIX64 " ", diff_ab[i - 1]);
        printf("\n");
*/
        base_type sum_c0_c2_carry = 0;
        base_type * sum_c0_c2 = diff_a; //Reciclamos memoria
        addc_blocks(sum_c0_c2,
                    c0,
                    c2,
                    sum_c0_c2_carry,
                    half_blocks * factor);
/*
        printf("sum_c0_c2 = ");
        for (unsigned int i = half_blocks * factor; i > 0; i--)
            printf("%.16" PRIX64 " ", sum_c0_c2[i - 1]);
        printf("\n");
*/
        base_type * c1 = diff_ab; //Reciclamos memoria.
        base_type carry = sign ? 1 : 0;
        if (!sign)
        addc_blocks(c1,
                    sum_c0_c2,
                    diff_ab,
                    carry,
                    half_blocks * factor);
        else
        subc_blocks(c1,
                    sum_c0_c2,
                    diff_ab,
                    carry,
                    half_blocks * factor);
        base_type c1_carry = carry ^ sum_c0_c2_carry ^ (sign ? 1 : 0);
        //delete [] diff_ab;
        delete [] sum_c0_c2;

        //resultado = c0_low  c0_high
        //                    c1_low   c1_high
        //                             c2_low   c2_high
/*
        printf("a0 = ");
        for (unsigned int i = half_blocks; i > 0; i--)
            printf("%.16" PRIX64 " ", a0[i - 1]);
        printf("\n");
        printf("a1 = ");
        for (unsigned int i = half_blocks; i > 0; i--)
            printf("%.16" PRIX64 " ", a1[i - 1]);
        printf("\n");
        printf("b0 = ");
        for (unsigned int i = half_blocks; i > 0; i--)
            printf("%.16" PRIX64 " ", b0[i - 1]);
        printf("\n");
        printf("b1 = ");
        for (unsigned int i = half_blocks; i > 0; i--)
            printf("%.16" PRIX64 " ", b1[i - 1]);
        printf("\n");

        printf("c0 = ");
        for (unsigned int i = half_blocks * 2; i > 0; i--)
            printf("%.16" PRIX64 " ", c0[i - 1]);
        printf("\n");
        printf("c1 = ");
        for (unsigned int i = half_blocks * factor; i > 0; i--)
            printf("%.16" PRIX64 " ", c1[i - 1]);
        printf("\n");
        printf("c2 = ");
        for (unsigned int i = half_blocks * factor; i > 0; i--)
            printf("%.16" PRIX64 " ", c2[i - 1]);
        printf("\n");
*/

        if (high)
        {
            //if (carry)   
            //    inc_blocks(c2 + half_blocks, carry,    half_blocks);
            if (c1_carry)
                inc_blocks(c2 + half_blocks, c1_carry, half_blocks);
        }

        carry = 0;
        addc_blocks<base_type>(r_data + half_blocks * 0,
                               r_data + half_blocks * 0,
                               c0,
                               carry,
                               min(half_blocks * 2,
                                   blocks_data - half_blocks * 0));
        if (carry &&
            blocks_data - half_blocks * 0 > half_blocks * 2)
            inc_blocks(r_data, carry, blocks_data - half_blocks * 0);

        if (blocks_data > half_blocks * 1)
        {
            carry = 0;
            addc_blocks<base_type>(r_data + half_blocks * 1,
                                   r_data + half_blocks * 1,
                                   c1,
                                   carry,
                                   min(half_blocks * factor,
                                       blocks_data - half_blocks * 1));
            if (carry &&
                blocks_data - half_blocks * 1 > half_blocks * 2)
                inc_blocks(r_data, carry, blocks_data - half_blocks * 1);
        }

        if (blocks_data > half_blocks * 2)
        {
            carry = 0;
            addc_blocks<base_type>(r_data + half_blocks * 2,
                                   r_data + half_blocks * 2,
                                   c2,
                                   carry,
                                   min(half_blocks * factor,
                                       blocks_data - half_blocks * 2));
            if (carry &&
                blocks_data - half_blocks * 2 > half_blocks * 2)
                inc_blocks(r_data, carry, blocks_data - half_blocks * 2);
        }


        delete [] c0;
        delete [] c1;
        delete [] c2;

        //printf("r_data = ");
        //for (unsigned int i = blocks_a + blocks_b; i > 0; i--)
        //    printf("%.16" PRIX64 " ", r_data[i - 1]);
        //printf("\n");
    }
    else
    {
        base_type h, l;
        l = mul<base_type>(a_data[0], b_data[0], h);
        base_type carry = zero;
        r_data[0] = addc<base_type>(r_data[0], l, carry);

        if (high)
        {
            r_data[1] = addc<base_type>(r_data[1], h, carry);

            if (carry && blocks_data > 2)
                inc_blocks<base_type>(r_data + 2, carry, blocks_data - 2);
                            
        }
    }
}
#else
///////////////////////////////////////////////////////////////////////////////

#error Esta version es mas rapida,
#error pero desgraciadamente tiene algun fallo en los calculos.

//PRE: blocks % 2 == 0
template<typename Base_Type>
void karatsuba(Base_Type * r_data,
               const unsigned int blocks_data,
               const Base_Type * a_data,
               const Base_Type * b_data,
               const unsigned int blocks)
{
    using base_type = Base_Type;
    const base_type zero = static_cast<base_type>(0);

    //auto max = [](unsigned int a, unsigned int b)
    //{return a > b ? a : b;};

    auto min = [](unsigned int a, unsigned int b)
    {return a < b ? a : b;};

    auto half = [](unsigned int n)
    {return n % 2 == 0 ? n / 2 : n / 2 + 1;};

    const bool high = blocks_data > blocks;
    const unsigned int factor = high ? 2 : 1;

    if (blocks > 1)
    {
        const unsigned int half_blocks = half(blocks);
        const base_type * a0 = a_data;
        const base_type * a1 = a_data + half_blocks;
        const base_type * b0 = b_data;
        const base_type * b1 = b_data + half_blocks;

        //Calculamos c0.
        base_type * c0 = new base_type[half_blocks * 2];
        zero_blocks<base_type>(c0, half_blocks * 2);
        requete_fma(c0, half_blocks * 2,
                    a0, half_blocks,
                    b0, half_blocks);

        //Escribimos c0 en el resultado.
        base_type carry = 0;
        addc_blocks<base_type>(r_data + half_blocks * 0,
                               r_data + half_blocks * 0,
                               c0,
                               carry,
                               min(half_blocks * 2,
                                   blocks_data - half_blocks * 0));
        if (carry &&
            blocks_data - half_blocks * 0 > half_blocks * 2)
            inc_blocks(r_data, carry, blocks_data - half_blocks * 0);

        //Calculamos c2.
        base_type * c2 = new base_type[half_blocks * factor];
        zero_blocks<base_type>(c2, half_blocks * factor);
        requete_fma(c2, half_blocks * factor,
                    a1, half_blocks,
                    b1, half_blocks);

        //Escribimos c2 en el resultado.
        if (blocks_data > half_blocks * 2)
        {
            carry = 0;
            addc_blocks<base_type>(r_data + half_blocks * 2,
                                   r_data + half_blocks * 2,
                                   c2,
                                   carry,
                                   min(half_blocks * factor,
                                       blocks_data - half_blocks * 2));
            if (carry &&
                blocks_data - half_blocks * 2 > half_blocks * 2)
                inc_blocks(r_data, carry, blocks_data - half_blocks * 2);
        }

        base_type sum_c0_c2_carry = 0;
        base_type * sum_c0_c2 = c2; //Reciclamos memoria
        addc_blocks(sum_c0_c2,
                    c0,
                    c2,
                    sum_c0_c2_carry,
                    half_blocks * factor);

        bool sign_a = false;// = a0 < a1;
        bool sign_b = false;// = b1 < b0;
        base_type * diff_a = c0; //Reciclamos memoria.
        base_type * diff_b = diff_a + half_blocks;
        base_type carry_a = 1;
        subc_blocks<base_type>(diff_a, a0, a1, carry_a, half_blocks);
        sign_a = carry_a ^ 0x1;

        base_type carry_b = 1;
        subc_blocks<base_type>(diff_b, b1, b0, carry_b, half_blocks);
        sign_b = carry_b ^ 0x1;

        const bool sign = sign_a ^ sign_b;

        if (sign_a)
        {
            carry_a = 1;
            neg_blocks<base_type>(diff_a, diff_a, carry_a, half_blocks);
        }
        if (sign_b)
        {
            carry_b = 1;
            neg_blocks<base_type>(diff_b, diff_b, carry_b, half_blocks);
        }

        base_type * diff_ab = new base_type[half_blocks * factor];
        zero_blocks<base_type>(diff_ab, half_blocks * factor);
        requete_fma(diff_ab, half_blocks * factor,
                    diff_a, half_blocks,
                    diff_b, half_blocks);
        delete [] diff_a;


        base_type * c1 = diff_ab; //Reciclamos memoria.
        carry = sign ? 1 : 0;
        if (!sign)
        addc_blocks(c1,
                    sum_c0_c2,
                    diff_ab,
                    carry,
                    half_blocks * factor);
        else
        subc_blocks(c1,
                    sum_c0_c2,
                    diff_ab,
                    carry,
                    half_blocks * factor);
        base_type c1_carry = carry ^ sum_c0_c2_carry ^ (sign ? 1 : 0);
        delete [] sum_c0_c2;


        if (c1_carry &&
            blocks_data > half_blocks * 3)
            inc_blocks(r_data + half_blocks * 3,
                       c1_carry,
                       blocks_data - half_blocks * 3);

        //Escribimos c1 en el resultado.
        if (blocks_data > half_blocks * 1)
        {
            carry = 0;
            addc_blocks<base_type>(r_data + half_blocks * 1,
                                   r_data + half_blocks * 1,
                                   c1,
                                   carry,
                                   min(half_blocks * factor,
                                       blocks_data - half_blocks * 1));
            if (carry &&
                blocks_data - half_blocks * 1 > half_blocks * 2)
                inc_blocks(r_data, carry, blocks_data - half_blocks * 1);
        }



        delete [] c1;
    }
    else
    {
        base_type h, l;
        l = mul<base_type>(a_data[0], b_data[0], h);
        base_type carry = zero;
        r_data[0] = addc<base_type>(r_data[0], l, carry);

        if (high)
        {
            r_data[1] = addc<base_type>(r_data[1], h, carry);

            if (carry && blocks_data > 2)
                inc_blocks<base_type>(r_data + 2, carry, blocks_data - 2);
                            
        }
    }
}
#endif
    
///////////////////////////////////////////////////////////////////////////////

//Bit mas significativo.
template<typename Base_Type>
inline
unsigned int msb_hint(const Base_Type * data,
                      const unsigned int first_search_block)
{
    using base_type = Base_Type;
    const base_type zero = static_cast<base_type>(0);
    const unsigned int bits_per_block = sizeof(base_type) * 8u;

    for (unsigned int i = first_search_block + 1; i > 0; i--)
    {
        if (data[i - 1] != zero)
        {
            base_type value = data[i - 1];

            unsigned int count = 0;
            base_type last = zero;
            while(true)
            {
                last = value & static_cast<base_type>(0xf);
                value >>= 4;
                if (value != zero)
                    count += 4;
                else
                    break;
            }
            static unsigned int msb[] =
                {0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4};
            return (i - 1) * bits_per_block + count + msb[last];
        }
    }

    return 0;
}

///////////////////////////////////////////////////////////////////////////////

//Bloque mas significativo.
template<typename Base_Type>
inline
unsigned int msblock(const Base_Type * data,
                     const unsigned int blocks)
{
    using base_type = Base_Type;
    const base_type zero = static_cast<base_type>(0);

    for (unsigned int i = blocks; i > 0; i--)
        if (data[i - 1] != zero)
            return i;

    return 0;
}

///////////////////////////////////////////////////////////////////////////////

//Nota: quotient y remainder deben ser inicializados a 0.
template<typename Base_Type>
void divide(const Base_Type * dividend,
            const unsigned int blocks_dividend,
            const Base_Type * divisor,
            const unsigned int blocks_divisor,
            Base_Type * quotient,
            const unsigned int blocks_quotient,
            Base_Type * remainder,
            const unsigned int blocks_remainder)
{
    using base_type = Base_Type;
    const base_type zero = static_cast<base_type>(0);
    const base_type one = static_cast<base_type>(1);
    const unsigned int bits_per_block = sizeof(base_type) * 8u;

    auto min = [](unsigned int a, unsigned int b)
    {return a < b ? a : b;};

    auto max = [](unsigned int a, unsigned int b)
    {return a > b ? a : b;};

    const unsigned int msblock_dividend = msblock(dividend,
                                                  blocks_dividend);
    const unsigned int msblock_divisor = msblock(divisor,
                                                 blocks_divisor);
    //Si el dividendo es 0...
    if (msblock_dividend == 0)
    {
        zero_blocks<base_type>(quotient, blocks_quotient);
        zero_blocks<base_type>(remainder, blocks_remainder);
        return;
    }

    base_type * rem;
    if (msblock_dividend > blocks_remainder)
    {
        rem = new base_type[msblock_dividend];
        copy_blocks<base_type>(rem, dividend, msblock_dividend);
    }
    else
    {
        rem = remainder;
        copy_blocks<base_type>(rem, dividend, msblock_dividend);
        zero_blocks<base_type>(rem + msblock_dividend,
                               blocks_remainder - msblock_dividend);
    }

    zero_blocks<base_type>(quotient, blocks_quotient);


    //Supondremos que el divisor no es 0.

    const unsigned int bits_per_iteration = bits_per_block / 2;

    //Obtenemos los primeros "bits_per_block / 2" bits del divisor.
    const signed int msb_divisor = msb_hint(divisor,
                                            msblock_divisor - 1);

    const unsigned int q = msb_divisor / bits_per_block;
    const unsigned int r = msb_divisor % bits_per_block;
    const unsigned int block_d = r == 0 ? q - 1 : q;
    const unsigned int shift = r == 0 ? 0 : bits_per_block - r;
    base_type hd = divisor[block_d] << shift;
    if (block_d > 0 && shift > 0)
        hd |= divisor[block_d - 1] >> (bits_per_block - shift);
    hd >>= bits_per_iteration;
    hd++;
    const base_type head_divisor = hd;


#if 0
const char * format = is_same<base_type, uint8_t>::value ?  "%.2"  PRIX8  " " :
                      is_same<base_type, uint16_t>::value ? "%.4"  PRIX16 " " :
                      is_same<base_type, uint32_t>::value ? "%.8"  PRIX32 " " :
                      is_same<base_type, uint64_t>::value ? "%.16" PRIX64 " " :
                      "%u ";
#endif
#if 0
printf("dividend = ");
for (unsigned int i = blocks_dividend; i > 0; i--)
   printf(format, dividend[i - 1]);
printf("\n");
printf("divisor = ");
for (unsigned int i = msblock_divisor; i > 0; i--)
   printf(format, divisor[i - 1]);
printf("\n");
printf("msb_divisor = %i\n", msb_divisor);
#endif

    unsigned int msb_first_search_block = msblock_dividend - 1;
#if 0
printf("msb_first... = %u\n", msb_first_search_block);
#endif

    while(msblock_divisor < msblock(rem, msb_first_search_block + 1) ||
          (msblock_divisor == msblock(rem, msb_first_search_block + 1) &&
          less_equal(divisor, rem, msblock_divisor)))
    {
#if 0
printf("rem     = ");
for (unsigned int i = msblock_dividend; i > 0; i--)
   printf(format, rem[i - 1]);
printf("\n");
#endif

        const signed int msb1 = msb_hint(rem, msb_first_search_block);
#if 0
printf("msb1 = %i\n", msb1);
#endif
        const unsigned int q = msb1 / bits_per_block;
        const unsigned int r = msb1 % bits_per_block;
        const unsigned int block = r == 0 ? q - 1 : q;
        const unsigned int shift = r == 0 ? 0 : bits_per_block - r;
        base_type h = rem[block] << shift;
        if (block > 0 && shift > 0)
            h |= rem[block - 1] >> (bits_per_block - shift);
        const base_type head = h;

        msb_first_search_block = block;

        const unsigned int q_diff_size = msb_first_search_block + 1;
        //const unsigned int q_diff_size = msblock_dividend;
        base_type q_diff[q_diff_size];
        zero_blocks(q_diff + 1, q_diff_size - 1);
        
        q_diff[0] = head / head_divisor;

        const signed int shift2 = msb1 -
                                  msb_divisor -
                                  static_cast<signed int>(bits_per_iteration);
#if 0
printf("shift2 = %i\n", shift2);
#endif

        if (shift2 <= 0)
        {
            q_diff[0] >>= -shift2;
            if (q_diff[0] == zero)
            {
                //Si seguimos ejecutando el bucle es porque
                //dividendo >= divisor.
                q_diff[0] = one;
            }
        }

#if 0
printf("q_diff  = ");
for (unsigned int i = q_diff_size; i > 0; i--)
   printf(format, q_diff[i - 1]);
printf("\n");
#endif

        const unsigned int p_size = msblock_dividend;
        base_type p[p_size];
        zero_blocks(p, p_size);
        fma2<base_type>(p,
                        p_size,
                        divisor,
                        msblock_divisor,
                        q_diff,
                        1);

#if 0
printf("p       = ");
for (unsigned int i = p_size; i > 0; i--)
   printf(format, p[i - 1]);
printf("\n");
#endif

        if (shift2 > 0)
        {
            base_type tmp[max(p_size, q_diff_size)];

            shift_left_big<base_type>(tmp, q_diff, shift2, q_diff_size);
            copy_blocks<base_type>(q_diff, tmp, q_diff_size);

            shift_left_big<base_type>(tmp, p, shift2, p_size);
            copy_blocks<base_type>(p, tmp, p_size);
        }

#if 0
printf("p       = ");
for (unsigned int i = p_size; i > 0; i--)
   printf(format, p[i - 1]);
printf("\n");
#endif

        //quotient += q_diff;
        unsigned int i = 0;
        while(q_diff[i] == zero) i++;
        base_type carry = zero;
        addc_blocks<base_type>(quotient + i,
                               quotient + i,
                               q_diff + i,
                               carry,
                               min(blocks_quotient,
                                   q_diff_size) - i);
        if (carry &&
            blocks_quotient > q_diff_size)
            inc_blocks<base_type>(quotient + q_diff_size,
                                  carry,
                                  blocks_quotient - q_diff_size);

        //rem -= p;
        i = 0;
        while(p[i] == zero) i++;
        carry = one;
        subc_blocks(rem + i,
                    rem + i,
                    p + i,
                    carry,
                    block + 1 - i);
    }
#if 0
printf("bloques divisor menores que bloques resto: %u\n",
       msblock_divisor < msblock(rem, msb_first_search_block + 1));
printf("bloques divisor iguales que bloques resto: %u\n",
       msblock_divisor == msblock(rem, msb_first_search_block + 1));
printf("divisor menor que resto restringido a %u bloques: %u\n",
       msb_first_search_block + 1,
       less_equal(divisor, rem, msblock_divisor));
#endif
    if (rem != remainder)
    {
        copy_blocks(remainder, rem, blocks_remainder);
        delete [] rem;
    }
#if 0
    base_type leches[blocks_dividend];
    memcpy(leches, remainder, blocks_dividend * sizeof(base_type));
    fma8<base_type>(leches, blocks_dividend,
                    quotient, blocks_quotient,
                    divisor, blocks_divisor);
    if (memcmp(leches, dividend, blocks_dividend * sizeof(base_type)) != 0)
        printf("mal dividido\n");
#endif

#if 0
printf("quotient = ");
for (unsigned int i = blocks_quotient; i > 0; i--)
   printf(format, quotient[i - 1]);
printf("\n");
printf("remainder = ");
for (unsigned int i = blocks_remainder; i > 0; i--)
   printf(format, remainder[i - 1]);
printf("\n\n");
#endif
}

///////////////////////////////////////////////////////////////////////////////
